Results (Julia)

using AlgebraOfGraphics
using AlgebraOfGraphics: density
using CairoMakie
using DataFrames
using DataFramesMeta
using Arrow
using PartialFunctions
using LaTeXStrings
using RCall
using Statistics
using beforerr


set_aog_theme!()

include("utils.jl")

:j0_k_norm => (log10 => L"Log Normalized Current Density ($J_A$)")

begin
    dir = "../data/05_reporting"

    j_events_low_fit = load()
    j_events_tau_20_fit = load(tau = 20)
    j_events_high_fit = load(ts = 0.12)

    j_events_low_der = load(method="derivative")
    j_events_high_der = load(ts = 0.12, method="derivative")

    w_events = load("$dir/events.Wind.fit.ts_0.09s_tau_60s.arrow")

    # add a label column to the dataframes
    j_events_low_fit[!, :label] .= "1 Hz (fit)"
    j_events_high_fit[!, :label] .= "8 Hz (fit)"
    j_events_tau_20_fit[!, :label] .= "1 Hz, 20s (fit)"
    j_events_low_der[!, :label] .= "1 Hz (derivative)"
    j_events_high_der[!, :label] .= "8 Hz (derivative)"

    # filter high time resolution events
    j_events_high_fit = @chain j_events_high_fit begin
        filter(:len => >(240), _)
    end

    j_events_der = vcat(j_events_low_der, j_events_high_der, cols=:intersect)

    # combine the dataframes
    j_events = reduce(vcat, [j_events_low_fit, j_events_high_fit, j_events_tau_20_fit])

    j_events = @chain j_events begin
        filter(:"fit.stat.rsquared" => >(0.95), _)
    end

    println("Number of events: ", size(j_events, 1))
end

┌ Warning: automatically converting Arrow.Timestamp with precision = NANOSECOND to `Dates.DateTime` which only supports millisecond precision; conversion may be lossy; to avoid converting, pass `Arrow.Table(source; convert=false)
└ @ Arrow /Users/zijin/.julia/packages/Arrow/Y6R1E/src/eltypes.jl:273

Number of events: 61879

begin
    log_axis = (yscale=log10, xscale=log10)

    j_limit = (10^-1, 10^1)
    j_norm_limit = (10^-2, 2)
end

(0.010000000000000002, 2)

Check the discontinuities properties with radial distance

datalimits_f = x -> quantile(x, [0.05, 0.97])

#16 (generic function with 1 method)

# NOTE: log axis for density is not working now
function plot_l_j_r(df)
    fig = Figure(size=(1000, 1000))

    plt = data(df) * mapping(row=:r, color=:label)
    plt = plt * density(datalimits=datalimits_f)

    plt1 = plt * mapping(l_map)
    plt2 = plt * mapping(l_norm_map)
    plt3 = plt * mapping(j_map)
    plt4 = plt * mapping(j_norm_map)

    axis = (; yscale = identity)
    
    # axis = (; yscale=log10, limits=(nothing, nothing, 10^-3, 1))
    # l_axis = (; yscale = log10, limits=(nothing, (10^-5, 10^-3)))
    # l_norm_axis = (; yscale = log10, limits=(nothing, (10^-3, 1)))

    grid1 = draw!(fig[1:5, 1], plt1, axis=axis)
    grid2 = draw!(fig[1:5, 2], plt2, axis=axis)
    grid3 = draw!(fig[1:5, 3], plt3, axis=axis)
    grid4 = draw!(fig[1:5, 4], plt4, axis=axis)
    legend!(fig[0, 1:end], grid1, titleposition=:left, orientation=:horizontal)

    fig
end

plot_l_j_r (generic function with 1 method)

plot_l_j_r(j_events)

function plot_l_j_r_hist(
    df;
    maps=(l_map, l_norm_map, j_map, j_norm_map),
    f_maps=(row=:r, color=:label),
    normalization=:pdf,
)
    fig = Figure(size=(1200, 1000))

    plt = data(df) * mapping(; f_maps...) * (visual(Lines) + visual(Scatter))
    plt *= histogram(datalimits=datalimits_f, normalization=normalization)

    plts = [plt * mapping(m) for m in maps]
    axis = (; yscale=log10)

    grids = [draw!(fig[1:5, i], p, axis=axis) for (i, p) in enumerate(plts)]
    legend!(fig[0, 1:end], grids[1], titleposition=:left, orientation=:horizontal)

    fig
end

plot_l_j_r_hist (generic function with 1 method)

plot_l_j_r_hist(j_events)

Figure 1

maps = (l_log_map, l_norm_log_map, j_log_map, j_norm_log_map)
plot_l_j_r_hist(j_events; maps = maps, normalization = :none)

function plot_l_j_r_hist(
    df;
    maps=(l_map, l_norm_map, j_map, j_norm_map),
    f_maps = (color=:r,),
    size = (1200, 800),
)
    fig = Figure(size=size)
    axis = (; yscale=log10)

    plt = data(df) * mapping(; f_maps...) * (visual(Lines) + visual(Scatter))
    plt_h = plt * histogram(datalimits=datalimits_f, normalization=:pdf)
    plts = [plt_h * mapping(m) for m in maps]
    grids = [draw!(fig[1, i], p, axis=axis) for (i, p) in enumerate(plts)]


    plt_h = plt * histogram(datalimits=datalimits_f)
    plts = [plt_h * mapping(m) for m in maps]
    grids = [draw!(fig[2, i], p, axis=axis) for (i, p) in enumerate(plts)]
    legend!(fig[0, 1:end], grids[1], titleposition=:left, orientation=:horizontal)

    fig
end

plot_l_j_r_hist(j_events_low_fit; maps=maps, f_maps=f_maps)

@rput j_events
@rput j_events_low_fit
@rput j_events_tau_20_fit
@rput j_events_der


R"""
library(ggplot2)
library(ggpubr)
library(patchwork)
source('utils.R')
"""

R"""
density_plot <- function(df,
  x,
  density.y = "density",
  facet.by = "r",
  color = "label",
  strip.position = "top")
{
  p <- ggdensity(
    df, 
    x = x, y = density.y,
    color = color,
    facet.by = facet.by,
    add = "median",
    alpha = 0
  ) %>%
  facet(facet.by, ncol = 1, strip.position = strip.position)
}

plot <- function(
  df,
  d1_ulim = 8000,
  d2_ulim = 40,
  d3_ulim = 8,
  d4_ulim = 1.2,
  xscale = "none",
  yscale = "log10",
  p1_xlim = c(10^2, 10^4),
  p1_ylim = NULL,
  p2_xlim = NULL,
  p2_ylim = NULL, # c(10^-3, 1)
  p3_xlim = NULL, p3_ylim = NULL,
  p4_xlim = NULL, p4_ylim = NULL,
  plot_func = density_plot,
  ...
) {

  p1 <- plot_func(
    filter(df, L_k < d1_ulim), "L_k", ...
  )
  
  p2 <- plot_func(
    filter(df, L_k_norm < d2_ulim), "L_k_norm", ...
  )
  
  p3 <- plot_func(
    filter(df, j0_k < d3_ulim), "j0_k", ...
  )
  
  p4 <- plot_func(
    filter(df, j0_k_norm < d4_ulim), "j0_k_norm", strip.position = "right", ...
  )

  p1 <- ggpar(p1,
      xscale=xscale, xlim = p1_xlim,
      yscale=yscale, ylim = p1_ylim
    ) + labs(x = lab_l)
  
  p2 <- ggpar(p2, 
      xscale=xscale, xlim = p2_xlim,
      yscale=yscale, ylim = p2_ylim, ylab=FALSE
    ) + labs(x = lab_l_norm)
  
  p3 <- ggpar(p3,
      xscale=xscale, xlim = p3_xlim,
      yscale=yscale, ylim = p3_ylim, ylab=FALSE
    ) + labs(x = lab_j)
  
  p4 <- ggpar(p4,
      xscale=xscale, xlim = p4_xlim,
      yscale=yscale, ylim = p4_ylim, ylab=FALSE
    ) + labs(x = lab_j_norm)
  
  tag_pool <- unique(df$r)
  
  p1 <- tag_facet(p1, open = "(a.", tag_pool = tag_pool)
  p2 <- tag_facet(p2, open = "(b.", tag_pool = tag_pool)
  p3 <- tag_facet(p3, open = "(c.", tag_pool = tag_pool)
  p4 <- tag_facet(p4, open = "(d.", tag_pool = tag_pool) + theme(
    strip.background = element_blank(),
    strip.text.x = element_blank()
  )

  p <- p1 + p2 + p3 + p4
  p <- p + plot_layout(nrow = 1, guides = "collect") & theme(legend.position='top')  
}
"""

R"""
# plot(
#     j_events_low_fit,
#     xscale = "log10",
#     d1_ulim = Inf, p1_xlim = c(1e2, 1e5), p1_ylim = c(1e-2, 1e0),
#     d2_ulim = Inf, p2_xlim = c(1e-1, 1e2), p2_ylim = c(1e-2, 1e0),
#     d3_ulim = Inf, p3_xlim = c(1e-2, 5e1), p3_ylim = c(1e-2, 1e0),
#     d4_ulim = Inf, p4_xlim = c(1e-3, 5e0), p4_ylim = c(1e-2, 1e0),
#     facet.by = NULL, color = "r"
# )

# save_plot("l_j_fit_xscale_log", width = 15, height = 5)

plot(
    j_events_low_fit,
    d1_ulim = Inf, p1_xlim = c(1e2, 1e5),
    d2_ulim = Inf, p2_xlim = c(1e-1, 1e2),
    d3_ulim = Inf, p3_xlim = c(1e-2, 5e1),
    d4_ulim = Inf, p4_xlim = c(1e-3, 5e0),
    facet.by = NULL, color = "r"
)

save_plot("l_j_fit", width = 15, height = 5)
"""

l_j_fit

l_j_fit_xscale_log

l_j_fit_xscale_log_count

Check the effect of time resolution and time window

R"""
# plot(j_events_der)
# save_plot("l_j_r_der", width = 15, height = 10)
# plot(
#     j_events,
#     p1_ylim = c(1e-5, 1e-3),
# )
# save_plot("l_j_r_fit", width = 15, height = 10)

plot(
    j_events,
    xscale = "log10",
    d1_ulim = Inf, p1_xlim = c(1e2, 1e5), p1_ylim = c(1e-2, 1e0),
    d2_ulim = Inf, p2_xlim = c(1e-1, 1e2), p2_ylim = c(1e-2, 1e0),
    d3_ulim = Inf, p3_xlim = c(1e-2, 5e1), p3_ylim = c(1e-2, 1e0),
    p4_xlim = c(1e-3, 5e0), p4_ylim = c(1e-2, 1e0),
)
save_plot("l_j_r_fit_xscale_log", width = 15, height = 10)
"""

l_j_r_fit

l_j_r_der

l_j_r_fit_xscale_log

Plot the median of discontinuity properties with radial distance

function q25(x)
    quantile(x, 0.25)
end

function q75(x)
    quantile(x, 0.75)
end

function stat_info(df)
    cols = [:L_k, :L_k_norm, :j0_k, :j0_k_norm]
    funcs = [median mean std q25 q75]
    group_cols = [:r, :label]
    
    @chain df begin
        groupby(group_cols)
        combine(cols .=> funcs)
    end
end

15-element Vector{Float64}:
  778.1492135059101
 1090.7343087931915
 1268.7318847139247
 1487.0954287365016
 1618.8611226298128
  734.1586231974484
  907.399617215176
 1195.9856252348263
 1357.2216902071673
 1378.3492860602692
  519.9710989502066
  697.3222159741033
  830.0951839086725
  999.2776112939654
 1098.9377386861229

arr= [1]
push!(arr, 2)
arr

true && 1

1

### Plot the median of discontinuity properties with radial distance
function plot_median_r(df; add_error=false)
    df_m = stat_info(df)

    fig = Figure()
    plt = data(df_m) * mapping(:r => r_lab, color=:label)

    mappings = [
        ("L_k", l_lab),
        ("L_k_norm", l_norm_lab),
        ("j0_k", j_lab),
        ("j0_k_norm", j_norm_lab)
    ]

    function plt_map(plt, col, label)
        plt *= mapping("$(col)_median" => label)
        lines_map = (visual(Lines) + visual(Scatter))
        errorbars_map = mapping("$(col)_q25", "$(col)_q75") * visual(Errorbars, alpha=0.1, whiskerwidth=10)

        maps = [lines_map]
        add_error && push!(maps, errorbars_map)
        
        plt * reduce(+, maps)
    end

    plts = [plt_map(plt, col, label) for (col, label) in mappings]
    popositions = [(1, 1), (1, 2), (2, 1), (2, 2)] # Define grid positions
    grids = [draw!(fig[pos...], plt) for (pos, plt) in zip(popositions, plts)]
    pretty_legend!(fig, grids[1])

    fig
end

plot_median_r(j_events)

plot_median_r(j_events; add_error=true)

Plot different radial distances in the same plot

function plot_l_r_one(df)
    fig = Figure(size=(1000, 500))

    plt = data(df) * mapping(col=:label, color=:r)
    
    grid1 = draw!( fig[1, 1:2], plt * mapping(:L_k_norm) * density(datalimits=((0, 50),)) )
    grid2 = draw!( fig[2, 1:2], plt * mapping(:L_k) * density(datalimits=((0, 5000),)) )
    legend!(fig[0, 1:end], grid1, titleposition=:left, orientation=:horizontal)

    fig
end

plot_l_r_one(j_events)

Figure 2

Plasma properties with radial distance

plt = data(j_events) * mapping(:radial_distance, :Alfven_speed, color=:label) * (visual(Scatter) + smooth())
plt |> draw

plt = data(j_events) * mapping(:radial_distance, jA_map, color=:label) * (visual(Scatter) + smooth())
plt |> draw

using Statistics

# groupby r and describe the data for each group 
# j_events |> @groupby(_.r) |> @map({r=key(_), j0_k=describe(_.j0_k), L_k=describe(_.L_k)})
@chain j_events begin
    groupby(:r)
    combine(:plasma_density =>  mean, :ion_inertial_length => mean, :b_mag => mean) 
end

5×4 DataFrame

Row	r	plasma_density_mean	ion_inertial_length_mean	b_mag_mean
	String	Float64	Float64	Float64
1	1.0	3.10512	152.217	3.75957
2	2.0	1.42267	242.722	2.10041
3	3.0	0.674452	317.706	1.45649
4	4.0	0.417696	434.665	1.3097
5	5.0	0.343857	569.222	1.27077

plt = data(j_events) * mapping(:radial_distance, :plasma_density, color=:label) * (visual(Scatter) + smooth())
plt |> draw

plt = data(j_events) * mapping(:radial_distance, di_map, color=:label) * (visual(Scatter) + smooth())
plt |> draw

Check discontinuities properties in relation to the local plasma properties

begin
    fig = Figure(size=(1000, 800))

    datalayer = mapping() * visual(Scatter, markersize=1, color=:white, alpha=0.1)

    begin
        layer = histogram(bins=range(1, 5, length=64), normalization=:pdf)
        plt = layer * mapping(di_log_map, l_log_map)
        plt1 = data(j_events) * mapping(layout=:r) * plt
        plt2 = data(w_events) * plt
    
        l_log_limit = ((1.5, 3.5), (1.5, 4.5))
        axis = (;limits=l_log_limit)
        draw!(fig[1:2,1:3], plt1, axis=axis)
        draw!(fig[1:2,4:5], plt2, axis=axis) 
    end

    # Current Density Panels
    begin
        layer = histogram(bins=range(-2, 3, length=64), normalization=:pdf)
        plt = layer * mapping(jA_log_map, j_log_map)
        plt3 = data(j_events) * mapping(layout=:r) * plt
        plt4 = data(w_events) * plt
    
        j_log_limit = ((-1, 3), (-2, 2))
        axis = (;limits=j_log_limit)
        draw!(fig[3:4,1:3], plt3, axis=axis)
        draw!(fig[3:4,4:5], plt4, axis=axis)
    end

    # Make ablines across different r panels
    # begin
        # ablines_data = (; intercepts=[-3,-1,1], slopes=[1,1,1]) 
        # lines = data(ablines_data) * mapping(:intercepts, :slopes) * visual(ABLines, linestyle=:dash)
        # draw!(fig[3:4,1:3], lines)
    # end

    Label(fig[0,1:3], "Juno in Different Radial Distances", fontsize=20)
    Label(fig[0,4:5], "Wind 11 Hz", fontsize=20)

    fig
end

Figure 3

begin
    fig = Figure(size=(1000, 800))

    pdflayer = density() * visual(Contour)
    # small scatter points
    datalayer = mapping() * visual(Scatter, markersize=3)

    # layer = pdflayer + datalayer
    layer = datalayer
    plt = layer * mapping(di_map, l_map)
    plt1 = data(j_events) * mapping(layout=:r) * plt
    plt2 = data(w_events) * plt

    l_limit = ((10^1, 10^5), (10^1, 10^5))
    axis = merge(log_axis, (;limits=l_limit))
    draw!(fig[1:2,1:3], plt1, axis=axis)
    draw!(fig[1:2,4:5], plt2, axis=axis)

    plt = layer * mapping(jA_map, j_map)
    plt3 = data(j_events) * mapping(layout=:r) * plt
    plt4 = data(w_events) * plt

    j_limit = ((10^-1, 10^3), (10^-2, 10^2))
    axis = merge(log_axis, (;limits=j_limit))
    draw!(fig[3:4,1:3], plt3, axis=axis)
    draw!(fig[3:4,4:5], plt4, axis=axis)


    Label(fig[0,1:3], "Juno in Different Radial Distances", fontsize=20)
    Label(fig[0,4:5], "Wind 11 Hz", fontsize=20)

    fig
end

Figure 4

begin
    fig = Figure(size=(1000, 800))
    datalimits_f = x -> quantile(x, [0.05, 0.95])

    begin
        layer = histogram(bins=24, datalimits=datalimits_f, normalization=:pdf)
        plt = layer * mapping(l_norm_log_map, j_norm_log_map)
        plt1 = data(j_events) * mapping(layout=:r) * plt
        plt2 = data(w_events) * plt

        axis = (;)
        draw!(fig[1:2, 1:3], plt1, axis=axis)
        draw!(fig[1:2, 4:5], plt2, axis=axis)
    end

    Label(fig[0, 1:3], "Juno in Different Radial Distances", fontsize=20)
    Label(fig[0, 4:5], "Wind 11 Hz", fontsize=20)

    fig
end

fig

Derivative method validation

fig = plot_l_j_r(j_events_der)
fig

Figure 5